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Abstract 

Howard Brenner"- •• • has recently proposed modifications to the Navier-Stokes equations that relate to a 
diffusion of fluid volume that would be significant for flows with high density gradients. In a previous pa- 
per/ we found these modifications gave good predictions of the viscous structure of shock waves in argon 
in the range Mach 1.0-12.0 (while conventional Navier-Stokes equations are known to fail above about 
Mach 2). However, some areas of concern with this model were a somewhat arbitrary choice of mod- 
elling coefficient, and potentially unphysical and unstable solutions. In this paper, we therefore present 
slightly different modifications to include molecule mass diffusion fully in the Navier-Stokes equations. 
These modifications are shown to be stable and produce physical solutions to the shock problem of a qual- 
ity broadly similar to those from the family of extended hydrodynamic models that includes the Burnett 
equations. The modifications primarily add a diffusion term to the mass conservation equation, so are at 
least as simple to solve as the Navier-Stokes equations; there are none of the numerical implementation 
problems of conventional extended hydrodynamics models, particularly in respect of boundary conditions. 
We recommend further investigation and testing on a number of different benchmark non-equilibrium flow 



1. Background 

A parameter which indicates the extent to which a local region of flowing gas is in thermodynamic equilibrium is the 
Knudsen number: 

Kn^j«- |Vp| , (1) 
L p 

where A is the mean free path of the gas molecules, L is a characteristic length of the flow system and p is a characteristic 
mass density. As Kn increases, e.g. for vehicles travelling at hypersonic speeds or at high altitudes, the departure 
of the gas from local thermodynamic equilibrium increases, and the notion of the gas as a continuum-equilibrium 
fluid becomes less valid. The Navier-Stokes equations (with standard no-slip boundary conditions) are, for example, 
typically confined to cases where Kn < 0.01. Their underlying constitutive laws for viscous stress tensor T and 
heat flux j^, i.e. Newton's law and Fourier's law respectively, may be derived from the Boltzmann equation using the 
classical Chapman-Enskog expansion in Kn to first order Extended, or modified, hydrodynamics models, such as the 
Burnett equations, are based on Chapman-Enskog expansions to higher orders in an attempt to extend the range of 
applicability of the continuum-equilibrium fluid model into the so-called 'intermediate-ZTn' (or 'transition-continuum') 
regime where 0.01 < Kn < 1. Extended expressions for T and include the same terms to first order in Kn contained 
in Newton's and Fourier's laws respectively, but with the addition of numerous, more complex terms that make them 
notoriously unstable and costly to solve. 

In 2004, ■ Howard Brenner proposed that the velocity u appearing in Newton's viscosity law is generally different 
from the mass velocity «„, appearing in the mass conservation equation: 

^ + V . (pu„,) = 0. (2) 

He subsequently related the two velocities by a diffusive volume fiux density j,, - u - u„ and proposed a constitutive 
model relating j,, to Vp, similar to Fick's law of mass diffusion. ' From this he derived modifications to the Navier- 
Stokes equations in which the governing transport equations of mass, momentum and energy remained unchanged but 
the constitutive equations for T and are augmented by additional terms. ^ The resulting equations have the form of an 
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extended hydrodynamics model, albeit one much less complex than the family of hydrodynamics models that includes 
Burnett, Grad, etc. 

The viscous structure of shock waves in gases provides an obvious test case for Brenner's modifications since 
they become increasingly significant as Vp increases. While it is accepted that the Navier-Stokes equations fail in nearly 
every respect in predicting correct shock structures above about Mach 2, they do reproduce the trends in experimental 
and molecular dynamics simulation data significantly better when Brenner's modifications are included, delivering an 
excellent match in the case of the inverse density thickness. ' It is of some concern, however, that shock solutions 
are nonphysical when the proportionality coefficient D^, used in the constitutive model for volume diff'usion, exceeds 
approximately the kinematic viscosity v - jj/p, where fi is the dynamic viscosity. Furthermore solutions were shown 
to be unstable when D,. > 1.45v. The imposed limit on D,. eff'ectively dictates the choice of D,, = v for which there is 
apparently no strong physical justification. 

Nevertheless, the results do partially support Brenner's original hypothesis — which is undoubtedly viewed with 
scepticism by some because it challenges the fundamental equations of fluid mechanics. Particular areas of criticism 
are that the underlying theory lacks a sound physical basis and is based on the notion of a diff'usive volume flux which 
is somewhat diflicult to conceptualise and for which constitutive models and their coefficients are untested. However, 
support can be found in the phenomenological GENERIC theory presented by Hans Christian Ottinger.' Ottinger 
questions why a diffusive transport term exists for both energy and momentum but not mass in the conventional Navier- 
Stokes equations and argues "something is missing", namely the ability of mass diffusion to produce entropy. When 
dissipative terms associated with mass density are identically zero, the GENERIC formulation arrives at the standard 
Navier-Stokes equations but, by including non-zero terms, a revised set of governing equations is derived that includes 
two velocities, similar to u,„ and u. The modifications are simply due to mass diffusion rather than the difficult concept 
of a diffusive volume flux. 

Brenner subsequently adopted^ the equations of Ottinger which differ from his original modifications particularly 
in that u appears not only in Newton's viscosity law but in the definition of momentum density itself. The purpose 
of this present paper is to provide additional argument in favour of the inclusion of mass diffusion and to examine its 
impact on the governing equations in detail. We assess the stability of the underlying equations and their ability to 
predict the structure of shock waves. 

2. Mass diffusion and conservation 

Thermal agitation causes molecules to travel from one region of a gas to another. Inequalities in molecular distribu- 
tion and thermal velocity tend to be smoothed by an inevitable net migration of molecules towards regions of lower 
molecular concentration and/or temperature. This is the process by which mass diffuses and, in ffie simplest case of a 
single specie gas, a diffusive flux jj occurs in the direction of negative density gradient, expressed through Pick's law 
simply as = ~-D,„'Vp where Dm is the coefficient of mass diff'usion (self-diffusion, in the case of a single specie gas). 
The conventional governing equations clearly omit the process of mass diffusion due to net migration of molecules by 
thermal agitation because they do not contain Pick's or any other constitutive model in the equation of conservation of 
mass. This is important in extended hydrodynamic modelling of hypersonic flows because the omission becomes more 
significant as Vp and, thus, Kn become larger 

Modelling of mass diffusion is, of course, a common feature in the analysis of multicomponent fluid systems. The 
usual approach is to retain the conservation equation for total fluid mass given by (|2|i and create additional conservation 
equations for the mass of individual gas species that each include a diffusive mass flux term.' However, for gas 
species, only - 1 equations of specie mass conservation are independent because the sum of all equations gives 
(|2]i. This means that mass diffusion is not modelled for one of the species, a statement that applies to the case of a single 
specie gas described in the previous paragraph. It happens because Q defines pu,„ to be the total mass flux density, 
i.e. the sum of the bulk, advective mass flux of the fluid and the net sum of diffusive mass flux of all consitutents of the 
fluid. In other words, m,„ constitutes a mean, or local mass average, velocity of all consitutents of the fluid' of which 
the advective velocity is only a part. 

The inability to model net mass diffusion and associated irreversible energy dissipation are not the only worrying 
consequences of combining advective and net diffusive fluxes into a single velocity «„,. It has been argued that, by 
doing this. Pick's law for a constituent is applied relative to the net flux of all constituents when it is only applicable 
relative to a frame of reference external to the fluid.' Also, velocity associated with mass diffusion has questionable 
physical significance because where the concentration of a given specie — > 0, its diffusion velocity — > oo.' 

Instead, let us spUt the total mass flux density pu„ - pu + where u is termed the advective velocity. If the 
diffusive flux is modelled by Pick's law, this leads to the following mass continuity equation: 

^+V.(pH)-V.(A„Vp)=0. (3) 
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This is the form of mass conservation equation proposed by Ottinger expressed in terms of u, not m,„. It is interesting 
to observe that the ratio R^c of diffusive mass flux to advective mass flux can be expressed as 



D,„Vp 



pu 



' (4) 



Aa yjySc Ma ' 

where: the Mach number of the flow Ma - \u\lc with c the speed of sound; the Schmidt number Sc - v/D„; y is the 
ratio of specific heats at constant pressure and volume; Kn is based on a MaxwelHan mean free path Xm - Ax yjyvic, 
with A,, = 16/(5 V2^) ~ 1.28. For argon gas, y = 5/3 and the coeflicient of self-diffusion of mass D,„ a; 1.32v,"' "' " 
so Sc - 0.76 and then R^c ~ Q.^Kn/Ma. In a planar shock in argon at Mach 4 upstream, Kn » 0.275 (see figure O 
and Ma » 1.0 local to the midpoint across the density profile, so R^c ~ 22%. The omission of mass diffusion from the 
governing equations will clearly lead to error at such a high R^/c- It could also be expected that R^c is high in regions 
of low speed and moderate density gradient such as boundary layers and wakes, both regions where the departure from 
non-equilibrium behaviour is most pronounced in hypersonic flows. ^ 

3. Momentum and energy conservation 

Brenner's original hypothesis was that Newton's viscosity law should be expressed in terms of u , not «„,. Along 
with his own supporting analytical and experimental evidence, there is also the argument that velocity in Newton's 
viscosity law represents deformation of fluid volume and therefore cannot be based on a mass average velocity m„,.^ 
The argument that velocity relating to mass diffusion has no physical significance' effectively precludes the use of 
in Newton's viscosity law. 

Ottinger additionally defines momentum density as pu, not pu,„, so that momentum relates purely to advective 
mass flux, not diffusive mass flux. This seems reasonable given that the advective flux is caused by mean translatory 
motion of molecules, associated with mechanical energy whereas the diffusive mass flux is caused by random motion of 
molecules associated with thermal energy. At the very least, if some momentum is attributed to mass diffusion, it must 
relate to a separate driving force independent of viscous forces associated with advective momentum. ' The resulting 
equation can be alternatively viewed as the governing equation for momentum in the absence of a (net) diffusive mass 
flux. Following these arguments, the momentum equation is as in the standard governing equations (ignoring body 
forces) but expressed in terms of «, not u,„: 

Du dipu) 

p-^^^+^-ipUmU)^V-P, (5) 

Dt at 

where the stress tensor P - T - pi (defined as positive in tension), p is pressure and / the unit tensor Newton's law is 
expressed asT -2p. dev(£)) + Kti{D)I where k is the bulk viscosity, the deformation gradient tensor D = symm(VM) = 
(1/2) I^Vm + (VM)^j and its deviatoric component dev(D) = D - (1/3) tr(£))/. Note that the material derivative DjDt is 
decomposed into the local rate of change djdt and the convective rate of change based on the local velocity of a fluid 
element that both advects and diffuses, i.e. «„,. 

The derivation of a conservation equation for energy within a continuum framework that includes mass diffusion 
is more challenging. The energy equation derived using GENERIC,' for example, contains an unconstrained phe- 
nomenological parameter that has to be determined by theory or simulation and confirmed by experiment. Derivation 
through physical argument alone requires careful accounting of contributions of mass diffusion to mechanical and ther- 
mal energies and associated work. As a first approximation, we equate the rate of change of total energy to mechanical 
and thermodynamic energy fluxes (ignoring internal sources): 

P^^^^ + '^-(pUmE)^V-(P-u)-V-j^, (6) 

Dt ot 

where the mechanical energy flux density iP-u) relates to the advective velocity u only and heat energy flux is attributed 
to conduction only by Fourier's law - -kVT where k is the thermal conductivity. The total energy per unit mass E 
represents all mechanical and thermal energy contributions. Initial simulations of the shock structure problem showed 
the temperature-density separation, discussed in section l631 was strongly underpredicted when the E included a kinetic 
energy due to advective mass flux only. Instead, much better predictions were obtained when the kinetic energy was 
due to the total mass flux such that E = e + |h„,P/2, suggesting that net mass diffusion contributes to an additional 
source of energy beyond the internal energy. 

4. Stability analysis 

The same stability analysis was undertaken on the set of governing equations proposed in this work that previously 
highlighted limitations of Brenner's original modifications.' Following the procedures described previously,' ''^ it is 
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first assumed that the gas is monatomic and calorically perfect with y - 5/3, Prandtl number fr = [yR/{y- l)](ij/k) - 
2/3, where R is the gas constant, and a- = 0. The governing equations from sections |2] and |3] are linearised in 1- 
dimension to produce the following non-dimensionalised perturbation equations: 



dcf> 



1 

1 1 
f 




(7) 



where 



1 dp' , 4 5m' , , 5dT' 

c - - — -z—, cr - --- — and q - --— — . (8) 

Sc ax' 3 ox' 2 ox' 

We assume a solution to (|7]i of the form 

= ^exp{i(wf' (9) 

where ^ is the amplitude of the wave, oj is its frequency and k its propagation constant. Equations (|7]i to Q can be 
combined to produce a set of linear algebraic equations of the form 

:R{u),k)4>^Q, (10) 

for which non-trivial solutions require 

det[^(w,A;)] = 0. (11) 
For our modified governing equations, (fTTT l yields the following characteristic equation: 

bios' + (23 + ()Sc-^)k^u^ - [lOk^ + (20 + 235c"')/t^]iw - [(15 + ASc-^)k'^ + 205c"'/t^] = 0. (12) 

If a disturbance in space is considered as an initial-value problem, k is real and o) - o),- + i<y, is complex. The form of 
(|9]l indicates that stability then requires w,- > 0. If a disturbance in time is considered as a problem of signalling from 
the boundary, w is real and k - kr + ikj is complex. For a wave travelling in the positive x direction, kr > 0, and stability 
then requires that kt < 0. For a wave travelling in the negative x direction, the converse is true: k,- < and stability 
requires A:, > 0. 

We examine temporal stability by solving (fT2l i numerically for oj for values of k in the range < A; < oo. 
Trajectories of oj are plotted in the complex plane in figure [Tl a). Stability was tested across a broad range of 0.2 < 
Sc < 1.0 and sets of trajectories are plotted at the two extremes. In both cases the trajectories all lie within the region 
CDj > 0, indicating stability for all k. 

We then turn to examine spatial stability by solving (T% numerically for k for values of ll) in the range < <y < oo. 
Trajectories of k are plotted in the complex plane in figure [TJb) for the same values of Sc as before. In both cases, the 
trajectories do not violate the stability condition. The results therefore show stable solutions for the modified Navier- 
Stokes equations presented in this paper 
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5. The shock structure problem and solution procedure 

The shock structure problem presented in this paper concerns the spatial variation in fluid flow properties across a 
stationary, planar, one-dimensional shock in argon. We define the flow as moving at an advective speed u in the 
positive x-direction, with the shock located at x = 0; the upstream conditions at x = -oo are super/hypersonic and 
denoted by a subscript '1', downstream conditions at x = +00 are denoted by a subscript '2'. Shocks were simulated 
for a range of upstream Mach number 1.2 < Mai < 1 1.0 with the problem initialisation and viscosity model detailed 
previously, ' but outlined briefly below. 

The simulations adopted the same thermodynamic models and coeflicients for argon used in earlier sections. The 
diffusive transport models were those previously described with the ratios of coefficients specified for argon, notably 
Pr = 2/3 and Sc - 0.76. The viscosity-temperature relation for argon was modelled by a power law of the form 
ju - AT'^ using an exponent s = 0.72 from independent experimental data. Since the results for this problem are 
historically presented in normalised form, the test problem was specified conveniently in a nondimensionalised form. 
The viscosity coefficient was set to A = 1, and in all simulations pi - Ti = 1 was specified at the upstream boundary. 
A gas constant R - y ' = 3/5 was chosen so that ci = 1 and, simply, ui - Mai for the particular simulation. At the 
downstream boundary, the normal gradient was specified as zero for all dependent variables except U2 whose value was 
specified using the Rankine-Hugoniot velocity relation to maintain the shock stationary and fixed within the domain. 

Simulations were performed using our solver, described in detail elsewhere, ' developed using the open source 
Field Operation and Manipulation (OpenFOAM) software.' A solution domain of 33Ami was used in all simulations, 
wide enough to contain the entire shock structure comfortably. Initial results were obtained using the conventional 
Navier-Stokes equations that converged on a mesh of 800 cells to within 1 % of the solution extrapolated to an infinitely 
small mesh size. The results presented in this paper were produced with a mesh of 2000 cells, corresponding to a 
mesh size of ~ 0.017 Ami- Numerical solutions were executed until they converged to steady-state, at which point the 
residuals of all equations had fallen 5 orders of magnitude from their initial level. 

Physical properties such as p and T vary continuously through the shocks from their upstream to their down- 
stream levels over a characteristic distance of a few mean free paths. Results presented in this paper are normalised 
between and 1 and denoted in the following by the superscript '*', against distance through the shock, nondimen- 
sionalised by Ami- Where possible results are compared with actual experiments""" rather than numerical Direct 
Simulation Monte Carlo (DSMC) data, since the latter requires certain assumptions relating to the form of the inter- 
molecular force law. 



6. Results 

Figure |2ta) shows the variation of p* and T* through a shock of Mach 2.84 calculated using the Navier-Stokes and 
modified Navier-Stokes equations. The experimental density profile of Torecki and Walenta' is also shown. It is clear 
that the shock layer predicted by the conventional Navier-Stokes equations is too thin, whereas the modified Navier- 
Stokes equations produce good agreement with the experimental data. The main region of disparity is upstream of 
the shock layer (left hand side in the figure) where the prediction trails out and is flatter than the experimental data. 




-6 -4 -2 2 4 6 ^-6 -4 -2 2 4 
Distance through shock [Ami ] Distance through shock [Ami ] 



Figure 2: Simulated and experimental profiles of a stationary shock at: (a) Mach 2.84; (b) Mach 9.0 
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Figure 3: Simulated and experimental inverse density thickness (Lp ') data, versus shock Mach number. 

Similarly, figure |2fb) shows the predicted profiles for a Mach 9.0 shock compared with experimental density data of 
Alsmeyer.^ Again, standard Navier-Stokes equations produce a shock profile which is too thin when compared with 
experiment. However, the modified Navier-Stokes equations produce excellent agreement with the experimental data. 

6.1 Inverse density thickness 

Apart from direct comparison of calculated and experimental shock profiles, there are other shock parameters for which 
experimental and/or independent numerical data is available. The principal parameter is the non-dimensional shock 
inverse density thickness, defined as: 

ip' = -^|VpU.. (13) 
Pi-P\ 

Comparing ([T]i and ( IT3T l it can be seen that, in the absence of a characteristic length scale L in an unconfined flow, the 
definition of Kn requires a characteristic dimension of a flow structure, in this case the actual thickness of the shock 
layer itself. Therefore ' has the interesting feature that it represents Kn for the shock structure case. Alsmeyer' re- 
ported the most comprehensive collection of experimental shock data, consisting of his own results and work published 
previously. Figure [3] shows L^' for argon shocks up to Mach 11, comparing simulation results and experimental data 
from Alsmeyer' and other sources.'-^ The Navier-Stokes equations predict shocks of approximately half the measured 
thicknesses over the entire Mach number range. As Lr^ indicates Kn, this poor agreement is expected because, over 
most of this Mach number range, Kn ~ 0.2-0.3, beyond the accepted Kn limit of application of the Navier-Stokes 
equations. However, results from the modified Navier-Stokes equations match well with experiment, suggesting that 
correct prediction of density gradient has been attained by correct modelling of mass diffusion. 

6.2 Density asymmetry quotient 

Agreement of predicted and experimental shock inverse density thicknesses is not the only measure of the success of 
a model. As L^' depends on the density gradient at the profile midpoint alone, it does not express anything about the 
overall shape of the profile. Instead, a second parameter that can be used to describe the shock profile, and for which 
experimental data is available, is the density asymmetry quotient Qp. This is a measure of how skewed the shock 
density profile is relative to its midpoint. It is defined for a 1 -dimensional profile of normalised density, p*, centred at 
p* - 0.5 on X = 0, as 

r p*(x)dx 

Qp = -j^ . (14) 

X [l-p*(x)]dx 

A symmetric shock would consequently have Qp - I, but real shock waves are not completely symmetrical about their 
midpoint. First, their general form is skewed a little towards the downstream. Then, the flattened, diffusive region, that 
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Figure 4: Asymmetry quotient (Qp) versus shock Mach Figure 5: Simulated and DSMC profiles of a stationary 
number. shock at Mach 1 1 

extends upstream of the shock profile, tends to increase Qp with increasing Mach number Figure|4]shows experimental 
data' in which Qp increases fairly linearly from ~ 0.9 at around Mach 1.5, through unity at around Mach 2.3, to ~ 1.15 
at Mach 9. 

Results from the Navier-Stokes equations do not agree well with experimental data: Qp > 1.0 at Mach 1 .2, and 
rapidly increases with Mach number with the profile sharpening downstream of the shock until by Mach 4 it levels off 
to 2p ~ 1 .4, compared to ~ 1 .03 from experiment. The modified Navier-Stokes equations similarly overpredict Qp at 
Mach 1 .2 but, with increasing Mach number, Qp quickly levels off at ~ 1.1 so that, by Mach 6, Qp matches well with 
experiment and the density profiles are very well predicted, e.g. for the profile at Mach 9 in figure|2b). 

6.3 Temperature-density separation 

In a shock, the density rises from its upstream value to its downstream value behind the temperature, due to the finite 
relaxation times for momentum and energy transport. Experimental data for this phenomenon is scarce due to the 
difficulty in measuring temperature profiles, but results from DSMC simulations provide such data. Figure |5] shows a 
comparison of profiles for a Mach 1 1 shock from our simulations and those calculated using DSMC.' DSMC clearly 
predicts a much larger separation distance between density and temperature profiles than conventional Navier-Stokes 
equations. The modifications to Navier-Stokes do increase the temperature-density separation, though not to the extent 
of the DSMC predictions. The temperature profile predicted by the modified Navier-Stokes equations is generally less 
diffusive than that predicted by DSMC, particularly at the upstream end. 

7. Conclusions 

It is accepted that the conventional Navier-Stokes equations fail to predict correct shock structures above about Mach 
2, where the flow falls within the intermediate-AT? regime. Brenner's modifications to Navier-Stokes improve the 
predictions of shock structures considerably ' but only with the somewhat arbitrary choice of diffusion coefficient 
D^, - V based on an upper limit above which the equations produce unphysical solutions and, at even higher D,., 
become unstable. 

Rather than basing the modifications to Navier-Stokes on the notion of a diffusive volume flux, we instead present 
modifications due to the inclusion of a diffusive mass flux. The resulting set of governing equations deliver a similar 
improvement over conventional Navier-Stokes in reproducing the trends in the experimental and DSMC data, and in 
the case of the inverse density thickness produce a very good match. The new model uses a known coeflicient for 
self-diffusion of mass for argon and the equation set does not exhibit the unphysical and unstable behaviour previously 
observed with Brenner's modifications. The new equation set is extremely easy to solve: it contains none of the higher 
order derivatives of the Burnett family of models, nor the second derivative of density contained within Brenner's 
original model; indeed, the addition of a diffusive term to the mass conservation equation makes the new set arguably 
easier to solve numerically than the conventional Navier-Stokes equations themselves. 
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While it is important not to draw strong conclusions based on just one test case, the results are generally en- 
couraging. Our future aims are: to test this model further on a number of benchmark cases ranging from high-speed 
flows encountered in hypersonics to specific studies of difliision phenomena such as thermophoresis, and to refine and 
develop the models accordingly. 

Acknowledgements 

We would like to thank Steve Daley of Dstl Farnborough (UK), Henry Weller of OpenCFD Ltd. (UK), Howard Brenner 
of MIT (USA) and Art Corey of Colorado State University (USA) for useful discussions. This work is funded in the 
UK by the Engineering and Physical Sciences Research Council under grants GR/T05028/01 and EP/D007488/1, and 
through a Phihp Leverhulme Prize for JMR from the Leverhulme Trust. 



8 



